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ABSTRACT 

Using radiation hydrodynamics simulations, we explore the evolution of circumplanetary disks 
around wide-orbit proto-gas giants. At large distances from the star (~ 100 AU), gravitational insta- 
bility followed by disk fragmentation can form low-mass substellar companions (massive gas giants 
and/or brown dwarfs) that are likely to host large disks. We examine the initial evolution of these 
subdisks and their role in regulating the growth of their substellar companions, as well as explore 
consequences of their interactions with circumstellar material. We find that subdisks that form in the 
context of GIs evolve quickly from a very massive state. Long-term accretion rates from the subdisk 
onto the proto-gas giant reach ~ 0.3 Jupiter masses per kyr. We also find consistency with previous 
simulations, demonstrating that subdisks are truncated at ~ 1/3 of the companion's Hill radius and 
are thick, with {h/r) of > 0.2. The thickness of subdisks draws to question the use of thin-disk approx- 
imations for understanding the behavior of subdisks, and the morphology of subdisks has implications 
for the formation and extent of satellite systems. These subdisks create heating events in otherwise 
cold regions of the circumstellar disk, and serve as planet formation beacons that can be detected by 
instruments such as ALMA. 

Subject headings: planetary systems: circumplanetary disks, planet formation; radiation hydrodynam- 
ics, ALMA 



1. INTRODUCTION 



At orbital radii > 100 AU, gravitational instability 
(GI) may be the principal formation mech anism for 
mass i ve planets a nd brown dwarfs (Stamatellos et al 



|2007[ |Boley"2009). GIs can occur in self-gravitating 
disks whenever the destabilizing effects of gravity be 
come comparable to the stabilizing effects of disk pres 
sure and shear (Toomre 1964). These instabili ties typ 
ically ma nifest themselves as spiral structure (Durisen 
et al [2007 ), which result in disk heating and mass redis- 
tribution through shocks and gravitational torques. In 
cases of very strong instability, spiral arms are unable 
to self-regulate and these arms can fragment into bound 
clumps. Because frag mentation is niost likely to occur at 
distances - 100 AU ([Rafikov" '20091 |Boleyl[2QQ9l IClarke 
|& Lodato||2009D, nascent fragments are a few AU in ra- 
dius and contain mostly molecular gas. Their large size 
makes these clumps susceptible to strong tidal effects, 
including complete destruction due to inward migration 
or due to pericenter p assages if on eccentric orbits (Bo- 
ley et al. 2010[ Nayakshin 2010). As fragments radiate 



away energy, they contract, leading to a rise in inter- 
nal temperature. Eventually, temperatures become high 
enough to cause H2 dissociation, which leads to a rapid 
collapse of the clump as energy goes into dissociation in- 
stead of increasing the gas pressure. At this stage, the 
low-angular momentum material will form a protoplanet 
or proto-brown dwarf that is no longer susceptible to 
tidal stripping by the host star. The high-angular mo- 
mentum material will be unable to fall directly on the 
core's surface, and an accretion disk will form. In this 
study, we focus on the initial evolution of a circumplan- 
etary/brown dwarf disk (forthwith subdisk) and its evo- 
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lution within a circumstellar disk in the context of an 
object formed by GIs. 

Characterizing the effects of a high-mass ratio binary 
(star/planet) on subdisk evolution is a necessary step to- 
ward understanding the growth of massive planets. Cir- 
cumplanetary disk evolution is also of interest for satellite 
formation, as subdisk structure, lifetimes, and size will 
impact the formation of large moons. A combination of 
analytics | Quillen & Trilhng 1998; Can up fc Ward||2002 



Wardfc Canup 2010; Martin &; Lubow 20 l l |) and simula- 
tions ([P^An^elo et al. 2002; Machida 2008| |Aylitte fc Bate 
20091 Martin fc Lubow 2011) have been used to explore 



the morphology and evolution of circumplanetary disks. 
General features that have been found include subdisk 
truncation at ~ 1/3 of the host planet's Hill radius and 
complex, three-dimensional gas flows. 

Most studies to date have focused on low-mass planets, 
and few of these stud ies use three-dimensional hy drody- 
namics (see Machida| [2QQ8l JAyli ffe fc Bate||2009f for ex- 



amples of simulations in 3D). Moreover, the inclusion 
of radiative transfer in these simulations is rare, which 
can play a s igniflcant role in regu l ating the evolution of 



the system ( [Ayliffe fc B"a:te||20Q9l 
we focus on the early evolution 01 



|2012D . In this study, 
newly- formed massive 
subdisks that are expected to form during massive planet 
and brown dwarf formation through disk fragmentation. 

This work uses three-dimensional radiation hydrody- 
namics simulations to explore the evolution of subdisks 
embedded in circumpstellar disks. The simulations are 
designed to address the following questions: (1) Is self- 
gravity likely to play a role in the evolution of subdisks 
born by disk fragmentation? (2) What is the role of the 
subdisks in regulating the mass growth of the proto-gas 
giant or proto-brown dwarf? (3) What are the observable 
signatures of nascent, wide-orbit proto-gas giants? 
This manuscript is organized as follows. In §2, we de- 
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scribe the computational method behind our radiation 
hydrodynamics simulations. The results are presented in 
§3, and their implications are discussed in §4 along with 
our conclusions. 

2. COMPUTATIONAL METHOD 

We aim to simulate and characterize the initial evo- 
lution of subdisks around newly-formed GI protoplan- 
ets. Simulations are run using CHYMERA, an Eulerian 
code that solves the equatio ns of hydro dynamics on a 



fixed, cylindrical grid (Boley et aL|[2QQ7a| ). The subdisks 
are evolved using r, (j)]z = 256, 512, 64 cells centered 
on the host planet. Each subdisk has a resolution of 
Ar,A(/),Az = 0.1 AU, 27rr/512, 0.1 AU. First, we in- 
vestigate an isolated subdisk (no circumstellar disk), re- 
ferred to as simulation (a), which is used for a base for 
comparisons. We then explore effects of a circumstellar 
disk on the evolution of the subdisk by allowing mass 
to flow onto the computational grid. This mass flux is 
varied with a gas volume density of 1 x 10~^^, 3 x 10~^^ 
and 1 X 10""*^^ gjcrv? for simulations (b), (c) a nd (d) , 
respectively. Further details are given in section 2.1.3 



Our simulations are initialized with a sub disk-to-planet 
mass ratio ~ 1 (3 Mj each). The proto-gas giant and cir- 
cumplanetary disk orbit the 1 Mq primary at 100 AU. 
The initial conditions (ICs) for the simulations are cre- 
ated using analytic approximations for the disk (axisym- 
metric, Keplerian rotation about the proto-gas giant, 
and non-self-gravitating vertical structure with a con- 
stant adiabatic index). Using these approximations the 
surface density is tune d to give an initially flat Toomre 
Q radial proflle fToomre |1964p with Q ~ 1.1. When 
Q < 1.7, di sks are suscept i ble to the formation of spiral 



arms (e.g., Durisen et al. 



2007)) and Q ~ 1 is a very 
unstable disE When the disk is run in complete isola- 
tion (no external potential from the host star), the disk 
is overall stable, although the approximations lead to 
vertical and radial oscillations. These out-of-equilibrium 
features are minor compared with the impact of includ- 
ing the host star's potential. For this reason, quiet er ICs 
are not developed, as the tidal potential (section 2.1.2) 
quickly alters the initial setup. 

Self-gravity, radiation transport, and the tidal poten- 
tial from the primary are included in our simulations . We 
use a modified version of the BDNL scheme from Boley| 
et al. (2007a) to handle the cooling of the gas. The algo- 
rithni invokes fiux-limited diffusion (FLD) in the radial 
and azimuthal directions, while ray tracing is used for 
the vertical direction where the majority of disk cooling 
occurs due to the large surface area of the photosphere. 
At very low optical depths, radiative transfer is in the 
free-streaming limit and the material cools slowly, as gas 
becomes inefficient at emitting and absorbing radiation. 
At very large optical depths, the high opacity also makes 
cooling inefficient. Ray tracing ensures that the tran- 
sition region between these two limits, which is where 
most of the disk cooling o ccurs, is followed ac curately. 



The equation of state from Boley et al. (2007b) is used, 
with the ortho-para hydrogen ratio fixed to 3:1. For the 
temperatures explored here, the gas is well approximated 
by an adiabatic index of 5/3. 

2.1. Modifications 



2.1.1. Radiative Transfer 

We use a subcycling method to control heating and 
cooling. When the radiative transport routine is called, 
we first calculate the divergence of the fiux for each cell. 
The ratio of the local internal energy density to the local 
divergence of the fiux gives a cell- wise radiative time step. 
We then compare the shortest radiative time step over 
all cells with the Courant time step. Ideally, the shortest 
radiative time step must be resolved in order to keep the 
algorithm stable. Time steps that are too large can pro- 
duce excessive heating or cooling of the gas, leading to 
unphysical results and numerical instability. We address 
this issue by allowing the radiative transfer to subcycle 
in cases where the radiative timescale is shorter than the 
hydrodynamic timescale, i.e., allow multiple calls to the 
radiative transfer routine for any one hydro step. We ap- 
ply the divergence of the fiux to the radiative time step 
in order to update the internal energy density for each 
subcycle. If the cooling/heating time has not reached 
the hydrodynamic timescale after 8 subcycles, we allow 
the gas to cool to no lower than the background tem- 
perature (see section 2.3). We also limit the heating for 
the last subcycle to a 10% increase in internal energy 
density. We calculate the optical depth once per cycle 
holding the value constant during the iterative cooling 
routine. This method provides numerical stability and 
reasonable speed while still attempting to resolve the ac- 
tual radiative time scale. 

As subdisks are not isolated, they will be affected by 
incident radiation from their environment. We include 
this effect by setting a minimum irradiation background 
temperature, as a function of distance from the proto- 
planet. This reffects the temperature required to balance 
the radiation received by an otherwise cold subdisk: 



V e ' ai 



D 



(1) 



Here, Rp is the host planet radius, which we set to 2 
D is the radial distance from the subdisk center, is the 
effective temperature of the host planet, which we take 
to be lOOOK, and Tb is the ambient background temper- 
ature due to the star and its environment, which we set 
to lOi^T. The distance dependence is based on a thin, 
irradiated disk. Tacc is the temperature corresponding to 
the luminosity that results from mass accretion from the 
subdisk onto th e prot oplanet, and is described in more 
detail in section [2.1.31 

2.1.2. Tidal and Indirect Potentials 

Simulations are performed in the frame of the planet, 
with the planet fixed to the center of the computational 
grid. The star's position is fixed in this reference frame at 
= x^e-x = —100 AU. To accommodate these choices, 
modifications are made to the potential field and the 
equations of motion. Let the total potential be 



(2) 



Here, the potential components from the star and planet 
are written as and <l>p, respectively. The indirect 
potential, ^ind^ captures the subdisk's infiuence on the 
planet, even though the planet is held fixed to the grid 
center (see |Michael fc Durisen,2010^ . The self-gravity of 
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the gas is given by ^sg^ and the rotational potential, 



R 



(3) 



takes into account the star-planet orbit. The vector x is 
the position of a grid cell relative to the planet. Be- 
cause the mass on the grid can be variable, we hold 
17^ = G{M^ + Mp)/xl fixed throughout the simulations. 
By including in the total potential, the equations of 
motion do not need to be changed to include the cen- 
trifugal terms. 

The Coriolis effect is not implicitly included in the po- 
tential terms, and must be added to the force calcula- 
tion. Let the radial momentum and angular momentum 
densities be represented by S and A, respectively. The 
Coriolis terms on our cylindrical grid are then 



S 
A- 



A - 2QvrrpAt 



(4) 
(5) 



for azimuthal and radial velocities and Vr^ respec- 
tively, gas mass volume density p, and cylindrical sub- 
disk radius r. With these updates applied during sourc- 
ing, the calculation can now be performed in the frame 
of the planet. 

2.1.3. Mass Flow Boundaries 

The disk boundary conditions are handled using both 
infiow and outfiow boundary cells. When material flows 
through the inner grid boundary (0.8 AU), mass is ac- 
creted by the protoplanet and included in the subsequent 
potential calculations. The process of accretion should 
generate an additional luminosity source that will irra- 
diate the subdisk. We account for this effect by mod- 
ifying our incident temperature profile, which sets the 
background radiation field at a given distance from the 
planet. This luminosity can be described by its temper- 
ature as 



GMpM 



(6) 



where a is the Stefan-Boltzmann constant. For each call 
to the radiative transfer routine, Tacc is calculated based 
on the current mass accretion from the subdisk M. This 
value is then used t o set t he background temperature as 
described in section [2.1.11 

The outer radial boundary uses outfiow conditions un- 
less the given boundary cell is set to fiux mass onto the 
computational grid. In the latter case, mass fiux is en- 
visaged to be part of the primary disk. One simulation 
(sim a) does not include any mass fiuxing, so all radial 
boundary cells use outfiow conditions. Simulations (b), 
(c), and (d) are run with incoming gas volume densities 
set to 1 X 10~^^, 3 X 10~^^, and 1 x 10~^^ g/cm^, respec- 
tively. For simplicity and overall stability, gas is fiuxed 
only up to half the vertical component of the domain. 
The disk's upper boundary is only an outfiow boundary. 
As will be shown below, the chosen volume densities lead 
to plausible surface densities for a fragmented protoplan- 
etary disk at 100 AU. 

3. RESULTS 
3.1. Mass Evolution 



We explore subdisk evolution by varying the amount 
of material that is allowed to fiow through the Hill re- 
gion in 4 different simulations (see Table fTl). The results 
show that the interactions between the subdisks and the 
circumstellar disks will have a large infiuence on the host 
planet. Figure 1 highlights this evolution for simulations 
(a) and (d), which shows the rate of mass fiowing from 
the subgrid onto the protoplanet (specifically the mass 
fiowing through 0.8 AU), the total mass on the compu- 
tational grid, and the cumulative mass that leaves the 
computational domain other than through accretion by 
the protoplanet. This mass lost profile shows an abrupt 
kink in simulation (d) and levels off soon after. This is a 
result of mass being fiuxed onto the grid, which quickly 
balances the mass that is lost. We find that all subdisks 
evolve away from their initial conditions within one orbit 
of the protoplanet about the star (1 kyr), which corre- 
sponds to about 7 orbits within the subdisk at 5 AU from 
the protoplanet (a little less than 1/2 the Hill radius). 
Most of the mass is accreted during this time period, 
reaching about 0.5 Mj for all simulations. After about 
1.5 kyr, the mass accretion rate settles down to a roughly 
steady state (Fig. [T]), with ~ 0.29 Mj/kyr in simulation 
(d), 0.14 in simulation (c), 0.11 in simulation (b), and 
0.10 in simulation (a). The subdisk in simulation (a) is 
not being continuously fed, so this rate can only persist 
for another ~ 6 kyrs. 

Simulation (d) retains the most massive subdisk, but 
even this disk quickly evolves away from its highly unsta- 
ble state. The minimum Q = Csn/ {-kGT^)^ for epicyclic 
frequency at the end of the simulation is ~ 1.7 for a 
narrow region at a subdisk radius ^2.5 AU. While this 
region is at the cusp of forming non-axisymmetric struc- 
ture due to GIs, structure appears to be dominated by 
the tidal potential, i.e., two distinct spiral arms aligned 
with the tidal distortion. 

3.2. Morphology 

Figures 3 and 4 show the end-state midplane tempera- 
ture and surface density images for each simulation. To 
define the limits of the disk, we also show in Figure 2 the 
cumulative mass profile as a function of distance from the 
protoplanet. In all simulations, the subdisk is truncated 
abruptly at ~ 4 AU, which is approximate 1/3 of the 
system's Hill radius, consistent with previous work. This 
truncation can be seen in the surface density images as 
well, and is highlighted by the circle in the midplane tem- 
perature images. The surface density images also show 
spiral arms in the subdisk resulting from the tidal poten- 
tial, as well as the spiral wakes, which penetrate from the 
circumstellar region all the way to the subdisk. Material 
enters the Hill sphere over a wide range of azimuths and 
creates shocks as circumstellar material meets the sub- 
disk. 

The interactions between the circumstellar disk and 
gas within the Hill region lead to shock structures that 
outline the wakes and circumstellar-subdisk boundary. 
These features are most pronounced in simulations (c) 
and (d). Figure 5 further highlights this, showing the 
velocity fiow of the gas around the Hill region in simula- 
tion (d). The heating of this gas could produce molecular 
walls as icy dust passes through the Hill sphere. While 
the subdisk itself also provides a region for heating ices, 
the molecular walls provide an additional signature of 
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an embedded clump. We highlight this in two different 
ways. First, we show in Figure 6 a radial cut at about 
6 o'clock for each panel in Figure 3. There is a rapid 
dropoff to the background temperature, followed by an 
isolated temperature spike in the Hill sphere shock re- 
gion, i.e., the potential molecular wall. While this region 
is still fairly cold, it will affect some ices such as CO. 
Moreover, we have only explored one location in the cir- 
cumstellar disk in these simulations, and similar features 
may be present at higher temperatures when the pro- 
toplanet is placed at smaller orbital radii for analogous 
simulatiuons. 

3.3. Prospects for ALMA 

In the previous section, we discussed how protoplan- 
ets/brown dwarfs in the outer regions of young, circum- 
planetary disks will alter the thermodynamics and the 
chemistry of the surrounding gas. Even if the protoplanet 
alone is difficult to detect, a relatively hot subdisk could 
significantly increase the continuum and line emission in 
localized regions of circumstellar disks. Morphological 
structures such as molecular walls, as defined above, can 
further reveal the presence of an embedded object. As 
shown in Figure 7, gas-phase CO column densities are 
~ 10^^ cm~^ in the molecular walls and ~ 10^^ cm~^ 
in the subdisk. Not only could the resulting emission be 
detected, but it could also be resolved. Molecular wall 
separations from the molecular subdisk are ~ few to 10 
AU wide, requiring at least O.''^! resolution. 

3.4. Vertical Structure: Thin vs Thick Circumplantery 

disks 

Figure 8 shows the end-state h/r for simulation (d), 
where the disk scale height h{r) = Cs{r)/Q{r). The 
sound speed and orbital frequency at a given r 
are determined using mass-weighted averages in the az- 
imuthal and vertical directions. For the simulated disk, 
the h/r 0.25 at the r - 1 AU and 0.5 at r - 5 AU. At 
even larger radii, h/r ~ 1, demonstrating that material 
flowing into the Hill sphere is not initially part of a disk, 
giving rise to complex three-dimensional structure. This 
is highlighted by a meridional slice through the substellar 
and antistellar locations (Fig. 9). 

These subdisks are thick (/i/r > 0.2), and may not be 
well- approximated by thin-disk theory {h/r < 0.1). The 
disparate concentrations of gas mass between subdisks 
and circumstellar disks have implications for understand- 
ing moon formation, as dust and ice coagulation, settling, 
and moon migration could be very different in the thick 
disk regime. These processes have largely been explored 
in the context of thin-disk theory (e.g., Canup & Ward 
2002), so simulations of gas-solid coupling and satellite 
migration in the context of disks with large scaleheights 
should be conducted. 

Why are these subdisks thick? Even though subdisk 
temperatures are low, subdisks are very hot compared 
with their orbital speed {h/r = Cs/[rQ]). At the outer 
boundary of the disk, the subdisk temperature would 
need to be 0.4 K for h/r - 0.1, weh below the CMB. At 
a subdisk radius of 1 AU, the disk temperature would 
need to be less than 10 K (below the stellar irradiation 
temperature) to have an /i/r ~ 0.1 for these simulations. 
At very small subdisk radii, the h/r may settle toward 



typical values in the thin-disk regime, as suggested by 
Figure 8. Once this regime is reached, then standard 
accretion disk theory may be applicable, but most of the 
subdisk mass will likely be contained in the thick-disk 
region. 

Finally, for reference, we calculate the effective a for 
a Shakura & Sunyaev (1973) a-disk that would be com- 
mensurate with the accretion rates seen in simulation 
(d). We calculate a using a steady-state solution, i.e., 
a ^ M /{Snc'^rt'^T^). For this estimate, we use the az- 
imuthal average for values at a subdisk radius ^ 1 AU, 
which yields a value of a = 0.0035. 

4. DISCUSSION AND CONCLUSIONS 

Planet growth and the formation of moon systems 
is dependent on the time evolution of circumplanetary 
disks. Radiation hydrodynamic simulations of these sys- 
tems in early stages of formation can give insight into the 
detailed physics of these processes. In this work, we have 
used radiation hydrodynamics simulations to investigate 
the initial evolution of circumplanetary disks under four 
different assumptions for mass inflow from a circumstel- 
lar disk. Despite the subdisks being initialized in an ex- 
tremely unstable state {Q ~ 1), the disks evolve rapidly 
away from strongly self- gravitating conditions. Simula- 
tion (d) does show a minimum Q ~ 1.7 at a subdisk 
radius 2.5 AU. This is at the boundary for the develop- 
ment of spiral distortions, but disk asymmetries appear 
to be driven by tidal features. While we cannot make 
a strong claim about the stability of the subdisk inte- 
rior to ~ 1 AU, the simulations presented here do argue 
against a sustained strongly gravitationally unstable sub- 
disk, where Q is near fragmentation conditions Q ^ 1.4 
(Mayer et al. 2004). 

After the initial phase of rapid evolution, each disk set- 
tles toward a rough steady state. Tides from the star pro- 
vide a strong enough perturbation on the subdisk to keep 
material flowing inward even though the subdisk is no 
longer significantly self-gravitating. The accretion rate 
of the subdisk onto the protoplanet is dependent on the 
mass of the circumstellar disk, although this dependency 
is shallower than linear for the regime explored here. In 
simulation (d), which is the most massive circumstellar 
disk used in these simulations, the mass accretion from 
the subdisk onto the planet is ~ 0.3 Mj per kyr. At this 
rate, we expect the protoplanet to double its initial mass 
in about 10 orbits. It should be noted that Gl-infected 
circumstellar d isks can evolv e significantly on timescales 
of 10^ yr (Boley et al. 2006), so some caution should be 
used when extrapolating the mass growth beyond this 
point. . We also note that an accretion mechanism will 
need to be available for the interior region of the subdisk 
(< 0.8 AU here) for subdisk mass to reach the proto- 
gas giant/brown dwarf, which is not captured in these 
simulations. 

Disk truncation is abrupt, and occurs at about 1/3 of 
the protoplanet 's Hill sphere. The specific angular mo- 
mentum distribution of material in the disk reaches a 
peak at about a subdisk radius of 3.5 AU, with a value 
^max ^ 6.5 X 10^^ cm^ s~^. A common parameterization 
for this value is often written as Jmax = i^i,R^^ where 
Rh ^ 12 AU for these simulations. The factor i is the 
angular momentum bias, and its value has been of consid- 
erable interest for creating circumplanetary disk models 
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(e.g., Ward and Canup 2010). Based on our simulation- 
derived values, ^ ~ 0.6. We note that can deviate 
considerably from the Keplerian value at the location of 
-^max, so Vt is taken directly from the simulation as done 
for the disk scale height calculation. 

Tidal effects from the host star on the subdisk, as well 
as material that falls into the Hill sphere from the circum- 
stellar disk, can produce gas temperature changes that 
create regions with gas-phase molecules that would oth- 
erwise be frozen on dust grains. These molecular walls, 
in addition to hotter regions of subdisks, could allow 
gas-phase processes in regions of circumstellar disks that 
would otherwise be too cold for certain interactions, as 
well as provide observational diagnostics. For example, 
we find CO column densities of ^ 10^^ to at least 10^^ 
cm~^, which could produce molecular emission observ- 
able with ALMA with potentially resolvable morpholog- 
ical structures. 

In addition to asymmetric structures, we find that 
subdisks have large h/r ratios, consistent with previous 
work that used very different methods. The tempera- 
ture required for these subdisks to be in the thin-disk 
regime is lower than realistic background temperatures, 
forcing subdisks to be thick. Three-dimensional simu- 
lations of subdisks that explore dust-solid evolution, in- 
cluding satellite-subdisk interactions, need to be explored 
in these thick disks. 

The answers to the three principal science questions of 
this work can be summarized as follows: (3) Subdisks 
create regions in outer disks that have higher tempera- 
tures than what would normally be present. These re- 
gions can lead large deviations in gas phase molecular 
abundances, which can be used for observational diag- 
nostics. Gas-phase chemistry that would normally be 
frozen out could be enhanced in these regions. 

We thank Paola Caselli, Peter Barns, and Eric Ford for 
discussions that improved this manuscript. We thank the 
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support was provided in part by the National Radio As- 
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the California Institute of Technology (Caltech) funded 
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Figure 1. Mass accretion history for the planet and subdisk for simulations (a) (no mass flux) shown as solid lines and simulation (d) 
(high mass flux) shown as dotted lines (see Table 1). The curves show the evolution of the subddisk mass (red), the mass lost from the 
computational grid (blue, isolated case only), the mass accreted by the planet (black), and the accretion rate onto the planet in Mj yr"-*^. 
We flnd that all subdisks evolve away from their initial conditions within one orbit of the protoplanet about the star (1 kyr), which 
corresponds to about 7 orbits within the subdisk at 5 AU from the protoplanet (about 1/2 the Hill radius). Most of the mass is accreted 
during this time period, reaching about 0.5 Mj for all simulations. After about 1.5 kyr, the mass accretion rate settles down to a roughly 
steady state, with ~ 0.29 Mj per kyr in simulation (d), 0.14 in simulation (c), 0.11 in simulation (b), and 0.10 in simulation (a). The 
subdisk in simulation (a) is not being continuously fed, so this rate can only persist for another ~ 6 kyrs. 



Table 1 

Initial Conditions 



All simulations 




Simulation 


Mass flux {g/cw?) 


Disk mass 


3 Mj 


(a) 





Host planet mass 


3 Mj 


(b) 


1 X 10-14 


Barycenter of subdisk orbit 


lOOAU 


(c) 


3 X lO-i'^ 


Disk radius 


lOAU 


(d) 


1 X 10-13 


Initial Toomre Q 


~1.1 






Host star mass 


IMo 







Circumplanetary Disks around Wide-Orbit Planets 




Figure 2. Cumulative mass radial profiles at 1.0, 2.0, and 3.0 kyrs in black, red and blue respectively, for simualations (a), (b), (c), and 
(d) (see Table 2). Vertical dotted lines represent the radial location of the Hill radius, 1/2 the Hill radius, and 1/3 the Hill radius. As the 
disk evolves, the largest fraction of mass remains within 5 AU. This is ~ 1/3 the Hill sphere of the host planet. 
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Figure 3. End-state midplane temperature contours (K) of the subdisk at 3 kyrs for simualations (a), (b), (c), and (d) (see Table 2). The 
disk-planet system's Hill radius is ~ 11.53 AU shown as black circles in each panel. The disk is truncated at Rhui/^- The regions of the 
disk where material is flowing onto the grid produce high temperature regions that could lead to ice desorption and rotational emission 
from gas phase molecules. 
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Figure 4. Surface density contour of embedded subdisk at 3 kyrs, in log g/cm?^ for simualations (a), (b), (c), and (d) (see Table 2). At 
radial distances just inside of 5AU, the material rapidly increases in density and velocity as the gravitational effects of the host planet start 
to dominate the flow of mass. The elongated shape of the subdisk is a result of tidal forces from the primary and dynamical interactions 
with circumplanetary material. 
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Figure 5. Surface density contour of embedded subdisk at 3 kyrs, in log g/cim?^ for simualations (d) (see Table 2). Shown in black are 
velocity vectors that trace the flow of mass into the subdisk from the surrounding circumstellar disk. 
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Figure 6. Radial cuts of the midplane temperature through the hot molecular wall at three different times, for simualations (a), (b), (c), 
and (d) (see Table 2). The red, green, and blue curves correspond to 1.0, 2.0, and 3.0 kyrs, respectively. Shock heating causes temperature 
variations along spiral arms and where material from the main disk material interacts with the subdisk. The inner disk is heated as mass is 
transferred inward by the spiral arms. The temperature variations in the subdisk and at the disk-subdisk boundary can cause desorption 
of some ices into the gas phase, increasing, e.g., CO abundances. Emission from these gas phase molecules may be observable with ALMA. 
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Figure 7. Gas-phase CO column densities (left) and the midplane CO number density in the gas phase relative to the total number of 
CO molecules (right). The gas-phase CO is largely confined to the molecular walls and in the warmer regions of the subdisk. Gas-phase 
CO can have very large column densities, potentially creating strong molecular line emission. 
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Figure 8. azimuthally averaged mass weighted scale height verses radial distance for simulation (d) at 3.0 kyrs. Subdisks are in the thick 
disk regime with aspect ratios {h/r) of > 0.2. 
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Figure 9. Volume density vertical cuts through the substellar and anti stellar regions of the sub disk. The disk aspect ratio shows that 
circumplanetary disks are in the thick disk regime, which will be important for future moon formation and migration studies. 



